#include "raster.hpp"
#include "rasterToField.hpp"
#include "rasterReader.hpp"
#include "meshTypeFunctions.hpp"
#include "meshTypes.hpp"
#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>

#include <sstream>

using Catch::Approx;
using namespace hms;
using namespace hms::raster;

TEST_CASE( "Creating raster" ){
	ArrayXXs values (2,3);
	values <<
		0, 2, 4,
		1, 3, 5;
	Array2s
		bottomLeftCorner {0,0},
		bottomLeftCenter {0.25, 0.25};
	
	Raster r1 { values.replicate(1,1), 0.5, bottomLeftCorner, true };
	Raster r2 { values.replicate(1,1), 0.5, bottomLeftCenter, false };
	REQUIRE( (r1.values() == values).all() );
	REQUIRE( (r2.values() == values).all() );
	REQUIRE( r1.bottomLeft().isApprox(bottomLeftCorner) );
	REQUIRE( r2.bottomLeft().isApprox(bottomLeftCorner) );
}

TEST_CASE( "Raster value at coordinates" ){
	ArrayXXs rasterValues (2,3);
	rasterValues <<
		0, 2, 4,
		1, 3, 5;
	
	Raster r { std::move(rasterValues), 2, {-2, -3} };

	REQUIRE( rasterValueAtCoords(r, Array2s{-3  ,-4}) == 0 );
	REQUIRE( rasterValueAtCoords(r, Array2s{-2  ,-4}) == 0 );
	REQUIRE( rasterValueAtCoords(r, Array2s{-1  ,-4}) == 0 );
	REQUIRE( rasterValueAtCoords(r, Array2s{-0.1,-4}) == 0 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 0  ,-4}) == 1 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 1  ,-4}) == 1 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 1.9,-4}) == 1 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 2  ,-4}) == 1 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 3  ,-4}) == 1 );

	REQUIRE( rasterValueAtCoords(r, Array2s{-3  ,-2}) == 0 );
	REQUIRE( rasterValueAtCoords(r, Array2s{-0.1,-2}) == 0 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 0  ,-2}) == 1 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 3  ,-2}) == 1 );

	REQUIRE( rasterValueAtCoords(r, Array2s{-3  ,-1.01}) == 0 );
	REQUIRE( rasterValueAtCoords(r, Array2s{-0.1,-1.01}) == 0 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 0  ,-1.01}) == 1 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 3  ,-1.01}) == 1 );

	REQUIRE( rasterValueAtCoords(r, Array2s{-3  ,-1}) == 2 );
	REQUIRE( rasterValueAtCoords(r, Array2s{-0.1,-1}) == 2 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 0  ,-1}) == 3 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 3  ,-1}) == 3 );

	REQUIRE( rasterValueAtCoords(r, Array2s{-3  ,0.99}) == 2 );
	REQUIRE( rasterValueAtCoords(r, Array2s{-0.1,0.99}) == 2 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 0  ,0.99}) == 3 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 3  ,0.99}) == 3 );

	REQUIRE( rasterValueAtCoords(r, Array2s{-3  ,1}) == 4 );
	REQUIRE( rasterValueAtCoords(r, Array2s{-0.1,1}) == 4 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 0  ,1}) == 5 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 3  ,1}) == 5 );

	REQUIRE( rasterValueAtCoords(r, Array2s{-3  ,10}) == 4 );
	REQUIRE( rasterValueAtCoords(r, Array2s{-0.1,10}) == 4 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 0  ,10}) == 5 );
	REQUIRE( rasterValueAtCoords(r, Array2s{ 3  ,10}) == 5 );
}

TEST_CASE("Raster to field values"){
	ArrayXXs rasterValues (2,3);
	rasterValues <<
		0, 2, 4,
		1, 3, 5;
	
	Raster r { std::move(rasterValues), 2, {-2, -3} };
	using MeshType = UniMesh;

	/* mesh origin is one unit upwards and to the right of raster origin.
	 * This way, ghost cells also get some original values. */
	MeshType mesh { uniMeshAs<MeshType>(3,4,1,1,false,-1,-2) };
	ArrayXXs fieldVals { rasterToFieldValues(r, mesh) };

	/* Raster with raster values:
	 *       |y
	 *       |
	 *  _____|_____
	 * |  4  |  5  |
	 * |_____|_____|
	 * |  2 _|__3__|____x
	 * |_____|_____|
	 * |  0  |  1  |
	 * |_____|_____|
	 * 
	 * mesh with cell indices:
	 *       |y
	 *       |
	 *     __|__ __ 
	 *  __|__|__|__|__
	 * |20|__|__|__|__|
	 * |15|__|__|__|__|_____x
	 * |10|__|__|__|__|
	 * |_5|__|__|__|__|
	 *    |_1|_2|_3|
	 * 
	 * overlay with raster values:
	 *       |y
	 *       |
	 *  __ __|__ __ 
	 * |__|_4|_5|_5|__
	 * |_4|_4|_5|_5|_5|
	 * |_2|_2|_3|_3|_3|_____x
	 * |_2|_2|_3|_3|_3|
	 * |_0|_0|_1|_1|_1|
	 * |__|_0|_1|_1|
	 * 
	 */

	// REQUIRE( fieldVals(0, 0) == 0 );
	REQUIRE( fieldVals(0, 1) == 0 );
	REQUIRE( fieldVals(0, 2) == 1 );
	REQUIRE( fieldVals(0, 3) == 1 );
	// REQUIRE( fieldVals(0, 4) == 1 );
	REQUIRE( fieldVals(0, 5) == 0 );
	REQUIRE( fieldVals(0, 6) == 0 );
	REQUIRE( fieldVals(0, 7) == 1 );
	REQUIRE( fieldVals(0, 8) == 1 );
	REQUIRE( fieldVals(0, 9) == 1 );
	REQUIRE( fieldVals(0,10) == 2 );
	REQUIRE( fieldVals(0,11) == 2 );
	REQUIRE( fieldVals(0,12) == 3 );
	REQUIRE( fieldVals(0,13) == 3 );
	REQUIRE( fieldVals(0,14) == 3 );
	REQUIRE( fieldVals(0,15) == 2 );
	REQUIRE( fieldVals(0,16) == 2 );
	REQUIRE( fieldVals(0,17) == 3 );
	REQUIRE( fieldVals(0,18) == 3 );
	REQUIRE( fieldVals(0,19) == 3 );
	REQUIRE( fieldVals(0,20) == 4 );
	REQUIRE( fieldVals(0,21) == 4 );
	REQUIRE( fieldVals(0,22) == 5 );
	REQUIRE( fieldVals(0,23) == 5 );
	REQUIRE( fieldVals(0,24) == 5 );
	// REQUIRE( fieldVals(0,25) == 4 );
	REQUIRE( fieldVals(0,26) == 4 );
	REQUIRE( fieldVals(0,27) == 5 );
	REQUIRE( fieldVals(0,28) == 5 );
	// REQUIRE( fieldVals(0,29) == 5 );
}

TEST_CASE( "Reading raster from istream" ){
	/* example from esri website, with adjusted values for ncols and nrows */
	std::stringstream pseudofile;
	pseudofile <<
		"ncols 4\n"
		"nrows 5\n"
		"xllcorner 378923\n"
		"yllcorner 4072345\n"
		"cellsize 30\n"
		"nodata_value -32768\n"
		"43 2 45 7 3 56 2 5 23 65 \n 34 6 32 54 57 34 2 2 54 6\n";
	
	RasterReader rr {
		{
			"nrows",
			"ncols",
			"xllcorner",
			"yllcorner",
			"xllcenter",
			"yllcenter",
			"cellsize",
			"nodata_value"
		}
	};
	rr.reverseY = false;

	Raster r { rr.read( pseudofile ) };

	ArrayXXs testValues(5,4);
	testValues.reshaped() = ArrayNs<20>{
		43, 2, 45, 7, 3, 56, 2, 5, 23, 65, 34, 6, 32, 54, 57, 34, 2, 2, 54, 6
	};

	REQUIRE( r.values().isApprox(testValues) );
	REQUIRE( r.nx() == 5 );
	REQUIRE( r.ny() == 4 );
	REQUIRE( r.cellsize() == Approx(30) );
	REQUIRE( r.bottomLeft().isApprox(Array2s{378923, 4072345}) );
	REQUIRE( r.noDataValue() == Approx(-32768) );
}